cpi uiaii 



Molecular kinetic analysis of a finite-time Carnot cycle 
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Abstract. - We study the efficiency at the maximal power rjmax of a finite-time Carnot cycle of 
a weakly interacting gas which we can regard as a nearly ideal gas. In several systems interacting 
with the hot and co ld reservoirs of the temperatures Th and T c , respectively, it is known that 
?7max = 1 — yJrT/Th which is often called the Curzon-Ahlborn (CA) efficiency -qcA- For the first 
time numerical experiments to verify the validity of t]ca are performed by means of molecular 
dynamics simulations and reveal that our 77 max does not always agree with ?7ca, but approaches 
t)ca in the limit of T c — > Th. Our molecular kinetic analysis explains the above facts theoretically 
by using only elementary arithmetic. 



Introduction. — Recently global warming has been 
a worldwide problem. Developing more efficient engines 
may help to solve such a problem. In physics, the effi- 
ciency of heat engines has been treated as a basic subject 
of thermodynamics. One of the most important results 
is the discovery of the Carnot efficiency which gives the 
upper limit of efficiency: vq = 1 — T c /T h , where T h and 
T c are the temperatures of the hot and cold heat reser- 
voirs, respectively. In spite of the high efficiency, rjc is 
usually realized only in the quasistatic limit. This means 
that the Carnot heat engine is useless as a real engine be- 
cause the power defined as output work per unit time is 
0. Real engines should work for a finite time and pro- 
duce a finite power. Therefore, the finite-time extension 
of the quasistatic heat engines is an important subject of 
thermodynamics. Curzon and Ahlborn [1,2] (see also [3]) 
considered such an extension of the Carnot cycle and de- 
rived a simple and beautiful result: the efficiency at the 
maximal power output is given by 



1 



(CA efficiency). 



(1) 



Several theoretical studies [4-9], ranging from the heat 
engine working in the linear response regime [6-8] to the 
heat engine working by a quantum mechanism [9] support 
the validity of Eq. (1). This implies that t]ca has some 
sort of universality independent of the model details. 

In spite of its importance, to our knowledge, no ex- 
periments have been carried out to verify the validity of 
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Eq. (1). Moreover, though in [1] the temperature differ- 
ences between the reservoirs and the working substance 
are taken as the parameters to maximize the power, they 
do not seem easily controllable. Thus, the CA efficiency 
Eq. (1) is, in our opinion, still controversial in these re- 
spects. 

In this paper, we consider a more natural extension of 
the quasistatic Carnot cycle as a model system by using 
a weakly interacting gas which we can regard as a nearly 
ideal gas. By means of molecular dynamics (MD) simula- 
tions, numerical experiments to verify the validity of the 
CA efficiency are performed for the first time. Our model 
also accepts theoretical analysis by using only elementary 
arithmetic. As shown later, we can reveal the validity and 
the limitation of the CA efficiency from that analysis. 

Model and simulations. — We consider the qua- 
sistatic Carnot cycle of an ideal gas first and then its 
finite-time extension. For simplicity, we here use the two- 
dimensional model. The usual quasistatic Carnot cycle 
of an ideal gas consists of four processes: (A): isothermal 
expansion process (Vi — ► Vz), (B): adiabatic expansion 
process (Vs — ► V3), (C): isothermal compression process 
(1/3 — ► V4), (D): adiabatic compression process (V4 — > Vi), 
where V^'s are the volumes of the cylinder at which we 
switch each of four processes (Fig. 1(a)). When we fix 
Th, T c , Vi and V2, we can easily determine the volumes V3 
and V4 since we assume an ideal gas as the working sub- 
stance. In fact, they are given by V3 = (T h /T c )V2 and 
V4 = (Th/T c )Vi for the two-dimensional case. In the case 
of a finite-time cycle, we assume that the right wall of the 
cylinder is a piston and moves back and forth at a constant 
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speed u. In our model, this u is taken as a unique param- 
eter to maximize the power, which is controllable unlike 
the parameters in [1]. We also assume that each process 
is switched at the same volume as in the quasistatic case. 

We have performed the two-dimensional event-driven 
MD simulations [10] as follows. We assume that N hard- 
disc particles with diameter d and mass m are confined 
into the two-dimensional cylinder with rectangular geom- 
etry and the collisions between hard-disc particles are per- 
fectly elastic. Defining (x,y) coordinates as in Fig. 1(b), 
we let the piston move along the x-axis at a finite con- 
stant speed u. Here, we express the x- length and the 
y- length of the cylinder as I and L, respectively. Then, 
the volume V, (i = l,---,4) of the cylinder at which 
we switch each of the four processes (Fig. 1(a)) becomes 
Vi = Lli, where k is the x-length of the cylinder at the 
switching volume Vi. If the process (A) begins at time 
t = 0, the volume V(t) of the cylinder at time t is given 
as V(t) = Ll(t) = L(ut + h) (0 < t < (l 3 - h)/u) in 
the expansion processes (A) and (B). V(t) in the com- 
pression processes (C) and (D) is also given as V(t) = 
L(-ut + 2l 3 - h) ((/ 3 - h)/u < t < 2(l 3 - h)/u). When 
a particle with the velocity v — (v x ,v y ) collides with 
the piston whose x-velocity is ±u, its velocity changes to 
v' = (—v x ± 2u,v y ). Therefore, the particle gives micro- 
scopic work m(|f| 2 — |i/| 2 )/2 = 2m(±uv x — u 2 ) against the 
piston. In the isothermal processes, to simulate the heat 
reservoirs, we set the thermalizing wall with the length 
S at the left bottom of the cylinder (see Fig. 1(b)). The 
thermalizing wall has the following feature [11,12]: When 
a particle collides with the thermalizing wall, its velocity 
stochastically changes to the value governed by the distri- 
bution function 




(— oo < v x < +oo,0 < v y < +oo, Ti (i = h in (A), c in 
(C))), where fee is Boltzmann constant. This thermaliz- 
ing wall may be understood as follows. Imagine a large 
particle reservoir thcrmalized at the temperature T; (i=h 
or c) instead of the thermalizing wall and assume that if 
a particle in the cylinder goes out into the particle reser- 
voir, another particle in the particle reservoir comes into 
the cylinder. From this consideration, we can see that the 
particles coming into the cylinder from the particle reser- 
voir obey the velocity distribution function proportional 
to the Boltzmann factor multiplied by v y . By normalizing, 
we can obtain the distribution function Eq.(2). As easily 
seen, this thermalizing wall guarantees that the particle 
velocities in the static system are governed by Maxwell- 
Boltzmann distribution with temperature Tj : 

The heat flowing from the thermalizing wall into the sys- 
tem can microscopically be calculated by the difference be- 



tween the kinetic energies before and after the collision on 
the thermalizing wall. We sum up the above microscopic 
heat during the simulation as well as the microscopic work. 
At the walls except the piston and the thermalizing wall, 
we adopt the reflecting boundary conditions for colliding 
particles. We have used N = 100 particles with d = 0.01 
and m = 1 in the system with L = 1, l\ = 1, I2 = 1.5, 
7h = 1, T c = 0.7, fee = 1 and S = 0.5. These parameters 
except T c are fixed in our all simulations and analysis be- 
low. As time progresses, thermodynamic variables should 
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Fig. 1: Schematic illustration of the model, (a) Pressure- 
volume (p-V) diagram of the quasistatic Carnot cycle for an 
ideal gas. (b) The piston moves at a finite constant speed u 
and the thermalizing wall with the temperature Ti (i=h, c) 
and the length S is set on the left bottom of the cylinder only 
in the isothermal processes. 

draw a steady cycle independent of initial states. Fig. 2 
shows the temperature- volume diagram for the steady cy- 
cle at u = 0.01 and u — 0.001, where k^T is determined 
as the kinetic energy per particle, assuming the principle 
of equipartition. From this figure, we can see that in the 
isothermal expansion (compression) process the temper- 
ature approaches a steady value lower (higher) than Xh 
(T c ) at u = 0.01. This result can easily be understood: 
If a heat engine is working at a finite u, heat should flow 
into the system at a finite rate to maintain the steady cy- 
cle. Therefore, the finite difference of the temperatures 



p-2 



Molecular kinetic analysis of a finite-time Carnot cycle 



u=0.01 
u=0.001 □ 
quasistatic 




Fig. 2: Temperature-volume (T-V) diagram for the steady cy- 
cle at it = 0.01 and u = 0.001. The data were obtained by 
averaging over 1000 cycles after transient 50 cycles in the MD 
simulations. The solid line is the quasistatic Carnot cycle of 
an ideal gas. We can see that as u becomes larger, the cycle 
obviously deviates from the quasistatic cycle. 



between the system and the heat reservoir is necessary. 
The cycle for u = 0.001 almost agrees with the quasistatic 
Carnot cycle of an ideal gas. This implies that our system 
of the hard-disc particles closely approximates an ideal gas 
system. 

We have also calculated the efficiency r\ = 
Wtotai/<3h,totai and the power P = Wtotai/V, where 
Wtotai is the total work against the piston, Qh, total is 
the total heat flowing into the system from the hot heat 
reservoir and r is the total time for the simulation. Fig. 3 
shows r] and P at various u. We have found that the 
maximal power is realized at u w 0.015. The corre- 
sponding efficiency 77 max (the efficiency at the maximal 
power) is about 0.18, which is close to the CA efficiency 
rycA = 0.163. 

Theoretical analysis. — To explain the above MD 
data, we construct the theoretical model using the ele- 
mentary molecular kinetic theory as below. We assume 
that even in a finite-time cycle, the gas relaxes to the 
uniform equilibrium state with a well-defined tempera- 
ture T very fast and the particle velocity v is governed 
by Maxwell-Boltzmann distribution /mb(», T). Wc would 
like to derive the time-evolution equation of T. The energy 
U of a two-dimensional equilibrium ideal gas is given by 
U = Nk B T. In the series of cycles, U can be changed by 
two factors: particle collisions with the thermalizing wall 
and the piston. Our strategy to derive the time-evolution 
equation of T is very simple: Counting the number of 
the particles colliding with the thermalizing wall and the 
piston and calculating the heat and the work from the dif- 
ference between the kinetic energies before and after the 
collisions. Firstly, we consider the effect of the thermal- 
izing wall. Since the number m of the particles with the 
velocity v (v y < 0) colliding with the thermalizing wall 
per unit time is given by iit = /mb(", T)S(— v y )N/V, the 
total number C of the particles colliding with the thermal- 
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Fig. 3: u dependence of the efficiency r\ and the power P. 
The bold dots are the values calculated by the MD simulations 
and the solid line is the values calculated by solving Eq. (9) 
numerically. The MD data were obtained by the simulations 
for 100-2000 cycles after 35-50 transient cycles. In the limit 
of u — * 0, the Carnot efficiency (r)c = 0.3) is realized and P 
vanishes. 



izing wall per unit time is calculated as 

SN \2vk B T 



C 



+ OC r-0 

dv x \ dvyTiT 



2-kV 



(4) 



The total energy of these colliding particles before the col- 
lisions is also given by 



/+oc />C 
dv x / 
-CO J — i 



m 2 3SNk B T 2nk B T 
dv y n T -v = \l-z-- (5) 



m 



Because the number of the reflecting particles is equal to 
the number of the colliding particles, the total energy of 
the particles after the collisions is calculated as 



C 



p-\-oo f-C 

/ dv v 

JO J 



dv x f(v,T^v 2 



iSNk B T t 2i:k B T 



4irV 



m 



(6) 



using Eq. (2). Therefore, the net energy transfer, namely 
the heat qi(t,T) flowing into the system per unit time in 
the isothermal processes (i = h in (A), c in (C)) is given 
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by 



3SNk B (Ti — T) 2nk B T 



4wV(t) 



m 



(7) 



Next, we derive the work against the piston by the col- 
liding particles in the expansion processes. To calculate 
the number of particles colliding with the piston, we con- 
sider the velocity distribution fMB{v x ,v y ,T) in the frame 
of the piston, where v x = v x — u and /mb(v x , v y , T) = 
/mb(£>x + u i v yi T) . The number tim of the particles with 
the velocity (v x ,v y ) colliding on the piston per unit time 
is n M = fyiB{v x ,Vy,T)Lv x N/V = /mb(v x , v y , T)v x N/l. 
Since a particle gives the work 2mu(v x — u) — 2muv x 
against the piston, the total work w e (t, T) against the pis- 
ton per unit time in the expansion processes becomes 



w e {t,T) 



/+oo /- + oo 

(IVy \ 
-oo Ja 



dv x 2muv x nM 



2muN \ A 2 T 3/2 Ff u 2 

l(t) \ 4 V 7T 2 



-I 



dv x [AVTv x 



(8) 



where A = y/2k B /m. The work w c (t, T) for the unit time 
in the compression processes is also obtained by changing 
u — > — u in Eq. (8). 

By the energy conservation law, the time evolution of T 
for each of four processes (A)-(D) is given by 



dT dT 
(A) : Nk B —=q h -w e , {B) : Nk B — = -w e 

(C) : Nk B ^- =q c - w c , (D) : Nk B ^- = -w c 



(9) 



Here, we have numerically solved the above Eq. (9) for 
the entire cycle. By using the final temperature of each 
process as the initial temperature of the next process re- 
peatedly, we can obtain the steady cycle of this heat en- 
gine. After reaching the steady cycle, we numerically cal- 
culate the efficiency i](u) = W(u)/Qh(u) and the power 
P{u) = W(u)u/ (2(^3 — li)), where Qh(u) is the heat trans- 
fer from the hot reservoir to the system, W(u) is the work 
output and 2(7 3 — h)/u is the time for one steady cycle. 
In Fig. 3, we plot the u dependence of 77 and P at the same 
parameters as in the MD simulations. From this figure, we 
can see that the correspondence between the MD data and 
the line calculated by solving Eq. (9) numerically is estab- 
lished qualitatively This implies that our assumption of 
fast relaxation to the equilibrium state is not so bad 1 . In 



1 Though one may see a discrepancy between the theory and the 
MD simulations in Fig. 3, we can see that a smaller S gives better 
agreement. This is because the speed giving the maximal power 
Wmax becomes small at a small S, which means that the gas is close 
to equilibrium and therefore meets our theoretical assumption that 
the gas always stays in the equilibrium state. This behavior of Umax 
will be confirmed in our analysis Eq. (17). 



0.8 



0.6 



0.4 



0.2 



numerical 

analytic • 

MD • 

CA efficiency 



0.2 



0.4 



0.6 



0.8 



Fig. 4: The efficiency at the maximal power J^max- The bold 
dots are the MD data, the solid line is the theoretical line 
calculated by solving Eq. (9) numerically, the dashed line is 
the the analytic expression Eq. (18) and the dotted line is the 
CA efficiency Eq. (1). 



Fig. 4, we compare the efficiency at the maximal power 
Vma,K = f?(umax), where u max is the speed giving the max- 
imal power, with the CA efficiency Eq. (1) at 7h = 1 and 
various T c . We have found that our rj max does not always 
agree with tjqa but tends to approach ?7ca as T c — ► Xh 
for both of the MD data and the numerical line. We have 
confirmed that this r; max behavior is common to the sys- 
tems with various parameters V\,V2,S, etc., though the 
data are not shown here. To explain this ?7 max behavior, 
we try to obtain the analytic form of ?7 max by solving the 
evolution equation of T in the following. 

As seen in Fig. 2, we can expect that T approaches 
a steady value Tjf in the isothermal expansion process 
(A). Then, 7^* is obtained as a solution of the equation 
dT/dt = in Eq. (9A). Because If = T h is realized in 
the quasistatic limit u — > 0, wc can expand Tjf by u as 



1 h 



•if u 2 



Eq. (9A), we can determine and 
up to 0(u 2 ) as 



0(u 3 ). Substituting Xf into 



and obtain Tjf 



rau 



t 4Lu r v^n , 
Th = Th -^s\— + k B 



35 J 



(10) 



If we assume that the relaxation to Tjf i s ver y f &s t> the 
heat flowing into the system during T(t) = is given by 



Q 



? = / 

Jo 



(h-h)/u 



g h (i,T h st ) dt 



Ql s - 2mNA^ (I + ^) uln ^,(1 ! ) 



using Eq. (7), where the quasistatic heat for ideal gas 
in the isothermal expansion process is defined as = 
Nk B T h ln(V 2 /Vi). Note that Q% -» Q£ s when we consider 
the quasistatic limit u — > 0. T c st and of the isother- 
mal compression process (C) can be obtained by replacing 
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T h , u, Vx and V 2 in Eqs. (10) and (11) with T c , -u, V 3 and 
V4, respectively. Firstly, we try to calculate 77 max by using 
and Qf above. By defining the work of one cycle W as 
W — Qh + Qf> we can calculate the efficiency 77 = W/Q^ 
and the power P = Wu/(2(l 3 — li)). The maximal power 
is realized at u = u max defined as a solution of dPj du = 0. 
Since u max is given by 



^max 



k B (T h -T c ) 



4m(i + 



27rfe E 



(12) 



T h + VT C [ 



Q^ and FF at u max are obtained as 
N 



Of(w) - ^bv^(v^+ V^)ln^, (13) 



Moreover, r; max = 77(u max ) is calculated as 



N V? 
-k B (T h -T c )ln^-. 



(14) 



(15) 



This is equal to the CA efficiency Eq. (1) though we ne- 
glect 0{u 2 ) in the calculation of <3 h * and Qf . Therefore, 
we may regard that this result gives a natural and mi- 
croscopic foundation of the original derivation of the CA 
efficiency Eq. (1). 

As seen in Fig. 4, however, the CA efficiency deviates 
from the MD data and the numerically calculated line. 
This is because there exists the heat transfer other than 
Qh and Q^, which may be missed in the original deriva- 
tion of Eq. (1) [1]. As seen in Fig. 2, the initial tempera- 
tures of isothermal processes are different from the steady 
values. This implies the existence of the additional heat 
transfer Q^ dd and Qc dd during the fast relaxation to the 
steady temperatures. Next, we repeat the similar deriva- 
tion of Eq. (15) by considering the effect of these addi- 
tional heat transfers Q^ dd and Q* dd . We define the total 
heat as Q h = Q^+Qf d and Q c = Qf+Qf d . Since we as- 
sume that the relaxation to the steady temperature is very 
fast, we can approximate the additional heat transfers as 
QJ dd = Nk B (T{? - T h ) and Qf d = Nk B (Tf - T c ), where 
Th and T c are the initial temperatures of the isothermal 
processes (A) and (C), respectively. If we assume that adi- 
abatic processes satisfy the relations T h = (V^/V^T^ and 
T c = (V^/V^Tjf which are the same as in the quasistatic 
case, we can obtain 



Ql d d = -Nk B 




(16) 



up to 0(u). Qc dd is given by changing Th <-> T c and 
u — > — u in Eq. (16). The work of one cycle W is defined 
as W = Qh' + Q? + <?h dd + Qf d - B y defining u max as a 
solution of dP/du = 0, we can obtain 



V2 8Lk B 
Vi + 3SAy/T£ 

T h \ I 1 , V 2 




From Fig. 4, we can see that Eq. (18) agrees with the 
MD data and the numerically calculated result very well 
due to the effect of the additional heat Qh dd and Q^ dd . 
To obtain the efficiency in the T c — ► Th limit, we set 
T c = T h - AT (AT < 1). Then, 7? max is given by 
77 max = AT/(2T h ) + C(AT 2 ) which is the same as the 
CA efficiency up to AT order. This result explains why 
our 77 max approaches tjca when T c — > Th. Very recently, 
similar ?7 max behavior has been observed also in the other 
types of the heat engines [13,14]. In the equilibrium limit 
of T c — > Th, the system may be regarded as being in the 
linear response regime. Therefore, our result is consistent 
with the CA efficiency proved by using the linear response 
theory [6]. 

Summary. — In this paper, we have studied the effi- 
ciency at the maximal power 7y max of a finite-time Carnot 
cycle of a weakly interacting gas which we can regard as a 
nearly ideal gas. Our model is a natural extension of the 
quasistatic Carnot cycle and has a piston moving back 
and forth at a constant speed u in the cylinder. We have 
used this u as a unique parameter to maximize the power. 
Since u is easily controllable, this model seems more natu- 
ral than the original Curzon-Ahlborn's model [1] . We have 
performed numerical experiments of this model by means 
of MD simulations to verify the validity of the Curzon- 
Ahlborn (CA) efficiency ?7ca for the first time and have 
found that our ?7 max does not always agree with tjca, but 
approaches tjca in the limit of T c — ► Th. Our molecular 
kinetic analysis can explain the above facts theoretically 
by using only elementary arithmetic. Especially, we have 
revealed that the difference between tjca and our ?7 max is 
due to the additional heat transfers which may be missed 
in the original derivation of tjca [1]- Though it is restricted 
in the equilibrium limit of T c — > Th, these results strongly 
support the validity of the CA efficiency from both of the 
experimental and theoretical points of view. We expect 
that our analysis in this paper will shed light on the mi- 
croscopic aspects of the finite-time extension of thermo- 
dynamics. 



We thank K. Nemoto and T. Nogawa for helpful discus- 
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